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Abstract 

The  multigrid  one  shot  method  for  optima]  control  problems,  governed  by  elliptic 
systems,  is  introduced  for  the  infinite  dimensional  control  space.  In  this  ca,sp  the  con¬ 
trol  variable  is  a  function  whose  discrete  representation  involves  increasing  number  of 
variables  with  grid  refiTiement.  The  minimization  algorithm  uses  Lagrange  muiiipliers 
to  calculate  sensitivity  gradients.  A  preconditioned  gradient  descent  algorithm  i.s  ac¬ 
celerated  by  a  set  of  coarse  grids.  It  optimizes  for  different  scales  in  the  representation 
of  the  control  variable  on  different  discretization  levels.  An  analysis  which  reduces 
the  problem  to  the  boundary  is  introduced.  It  is  used  to  approximate  the  two  level 
asymptotic  convergence  rate,  to  determine  the  amplitude  of  the  minimization  stej), 
and  the  choice  of  a  high  pass  filter  to  be  used  when  necessary.  The  effectiveness  of  the 
method  is  demonstrated  on  a  series  of  test  problems.  The  new  method  enables  the 
solutions  of  optimal  control  problems  at  the  same  cost  of  solving  the  corresponding 
analysis  problems  just  a  few  times. 
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1  Introduction 


Numerous  computational  methods  have  been  developed  for  predicting  the  performance  of 
physical  systems.  For  engineering  design  purposes  a  modification  of  the  system  configura¬ 
tion  that  results  in  optimal  performance  is  required.  However,  computations  of  large  scale 
optimal  control  problems  are  extremely  time  consuming  and  in  many  cases  not  practical. 
The  effort  to  overcome  such  computational  difficulties  is  done  in  the  direction  of  devel¬ 
oping  faster  computers,  on  the  one  hand,  and  in  improving  the  performance  of  existing 
algorithms,  on  the  other  hand. 

An  important  and  difficult  class  of  optimal  control  problems  are  optimal  shape  design 
(OSD)  of  aerodynamic  systems  [2,  8,  9,  10,  11|;  for  example,  the  design  of  a  wing  shape. 
Under  certain  assumptions,  OSD  problems  can  be  reduced  to  simpler  optimal  boundary 
control  (OBC)  problems  using  the  small  disturbance  approximation  [6].  The  resulting 
problems  involve  a  fixed  physical  domain  with  control  variables  that  are  defined  as  boundary 
data.  The  problem  is  to  minimize  a  cost  function  under  certain  constrains,  which  are  a  set 
of  PDFs  called  “state  equations”.  The  cost  function  is  defined  to  measure  the  performance 
of  the  physical  system. 

A  standard  solution  process  involves  an  iterative  algorithm  where  each  iteration  is  com¬ 
posed  of  two  steps.  First,  the  control  variables  are  updated  with  the  “sensitivity  gradients” 
which  are  the  gradients  of  the  cost  function  with  respect  to  the  control  variables.  Then  the 
state  variables  are  updated  by  solving  the  constraint  equation  with  the  new  values  of  the 
control  variables.  The  repeated  solution  of  the  constraint  PDE  makes  this  computation 
extremely  time  consuming  and  in  some  cases  not  practical. 

Several  methods  were  developed  to  calculate  the  sensitivity  gradients.  Among  them 
is  the  “adjoint  method”  [5,  9,  11].  In  the  adjoint  formulation,  a  Lagrangian  is  defined 
together  with  Lagrange  multipliers,  which  are  also  called  “costate  variables”.  Costate  and 
design  equations  are  derived  from  the  variation  of  the  Lagrangian,  and  together  with  the 
state  equation  form  the  necessary  conditions  for  a  minimum.  The  sensitivity  gradients  are 
the  design  equation  residuals  calculated  with  the  solutions  of  the  state  and  costate  PDFs. 
The  adjoint  method  was  first  applied  to  aerodynamic  design  by  A.  .Jameson  in  1988  [9].  A 
multigrid  (MG)  solver  was  used  to  accelerate  the  convergence  of  the  solution  of  the  state 
and  costate  equations.  This  reduces  the  computational  cost  of  each  optimization  step  to 
0{N)  operations  (where  N  is  the  number  of  computational  grid  points)  but  does  not  reduce 
the  number  of  iterations  to  reach  the  minimum. 

Originally  MG  methods  were,  developed  to  accelerate  the  convergence  rate  of  the  nu-  — 
merical  solutions  of  PDFs  [3,  7,  14].  A.  Brandt  suggested  in  1984  [4,  page  1 19]  to  apply  MG 
methods  for  optimization  problems  in  the  framework  of  a  full  multigrid  (FMG)  algorithm 
where  the  optimization  problem  should  be  solved  on  coarse  levels  and  interpolated  to  finer 
levels  until  the  finest  level  is  reached.  It  is  further  suggested  in  [4]  to  treat  the  optimization 
problem  on  all  levels  where  on  the  finer  grids  the  optimization  step  should  be  done  locally  — 
if  possible  and  the  smooth  corrections  of  the  error  should  be  done  during  the  coarse  grid  __ 
correction.  The  whole  problem  should  be  solved  in  one  application  of  the  FMG  solver.  ' 


□ 
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The  adjoint  and  M(J  methods  were  romhined  to  solve  an  optimal  control  proldem,  in 
“one  shot",  for  the  finite  dimensional  control  by  S.  Ta’asan  in  1991  [12].  In  ihe  finite  dimen¬ 
sional  approach,  the  control  variables  are  repres«*nted  as  a  finite  sum  of  some  preassip,ned 
base  functions.  The  main  idea  in  [12]  is  to  represent  the  state,  the  costate  and  the  design 
equations  on  coarser  grids  with  the  full  approximation  scheme  (FAS )[;!].  It  is  shown  in  [12] 
that  in  general  the  use  of  Lagrange  multipliers  is  essential  to  achieve  acceleration  of  tin*  fine 
grid  solution  process  by  coarser  grids.  The  algorithm  optimizes  the  control  varialdes  on 
coarse  grids,  and  thus,  eliminates  the  repeated  .solution  of  fine  grid  e(piati<jn.s  in  every  op¬ 
timization  step.  The  one  shot  algorithm  was  applied  successfully  to  the  small  disturbance 
approximation  of  an  aerodynamic  wing  design  prolilem  in  a  subsonic  flow  by  S.  la'asan, 
G.  Kuruvila  and  M.  D.  Salas  (1992)  [Id].  The  performance  of  the  algorithm  in  [LI]  was  a 
substantia!  improvement  in  terms  of  computational  cost.  However,  the  perfoiinance  of  the 
finite  dimensional  one  shot  algorithm  depends  on  on  the  choice  of  Irase  functions  and  on 
the  level  on  which  the  different  control  variables  were  optimized. 

In  this  paper  we  extend  the  multigrid  one  shot  algorithm  to  the  infinite  dimensional 
control.  We  introduce  an  analysis  which  reduces  the  problem  to  the  boundary.  The  analysis 
is  used  to  to  determine  a  minimization  step  which  reduces  mainly  the  high  frequency  errors 
in  the  control  variables.  In  elliptic  systems  such  a  minimization  step  requires  an  uj)date  of 
the  state  and  costate  solutions  only  in  a  local  area  neighboring  the  bounilary.  Ba.sed  on 
the  above,  two  level  analysis  is  done  to  approximate  the  convergence  performance  of  the 
algoritbtn  for  a  given  problem  and  discretization.  Computational  demonstrations  of  the 
algorithm  are  given  for  a  set  of  test  problems  in  which  the  PDE  constraint  is  elliptic. 

2  A  Single  Grid  Algorithm  for  the  Solution  of  Opti¬ 
mal  control  Problems 

2.1  Problem  Definition 

Let  n  be  a  bounded  open  .set  of  with  smooth  boundary  F  and  let  <p  be  a  real  valued 
function  on  Q.  Let  1/  and  W  be  Hilbert  spaces  of  real  valued  functions  which  are  defined 
on  r  and  Q  respectively. 

The  problem  is  to  find  the  “control  variable”,  u  £  U,  and  the  “state  variable”,  d)  £  W. 
such  that  a  given  cost  function,  F{u,4>{u)),  defined  on  U  x  W,  will  be  minimized.  Here 
satisfies  an  elliptic  PDE  which  is  defined  on  D  and  will  be  referred  to  as  the  “analysis 
problem”  or  the  “state  equation”: 

miiiuei/ /'"(«,  ^(m))  07?.  r  .  . 

L{<f),u)  =  0  on  Q 

Note  that  the  control  variable  is  defined  on  tlie  boundary  F,  therefore  an  “optimal  boundary 
control”  (OBC)  problem  is  considered. 
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2.2  Derivation  of  the  Necessary  Conditions  for  a  Minimum 

We  apply  the  adjoint  method  to  the  optimal  boundary  control  problem  (2.1  j.  The  varialile 
space  is  enlarged  by  adding  Lagrange  multiplier  functions  or  costate  varialjles  denoted  by 
A.  A  Lagrangian  is  defined  to  be  the  sum  of  the  functional  and  a  linear  term  in  the  custat(> 
variables  which  vanishes  as  the  constraint  equation  is  satisfied; 

E{<i>,X,u)  =  F{u,4>)  -  {X,L{(p,  u)) .  (2.2) 

A  perturbation  of  the  Lagrangian  with  respect  to  all  the  variables  independently,  i.e..  state, 
costate  and  control,  results  in  a  variation  of  the  Lagrangian: 

6  *—  4>  -{■  £<t> 

Ae-A  +  eA  (2.;3) 

u  e-  u  -f  eii 

with  A  6  ii  ^  U  and  e  is  a  small  real  parameter.  The  variation  of  the  Lagrange 

function,  bE,  in  the  first  order  approximation  in  e,  is  given  in  the  following  form; 

8E  =  -e  u)X  +  F^)  -  £  (a,  L{4>.  u))  +  £  (u.  L;(<?>,  u)A  +  F,)  (2.4) 

where  LJ,  and  L*  are  the  adjoint  operators  of  and  respectively.  The  requirement 
that  the  first  approximation  terms  vanish  results  in  the  necessary  condition  for  a  minimum 
which  will  be  referred  as  the  state,  the  costate,  and  the  design  equations: 

State :  = 

Costate  :  Z/J(0,  u)A  +  F^(^^»,  u)  =  0  (2.5) 

control  :  L*(0,  u)A  +  F„(^,  u)  =  0. 

From  here  on  we  will  use  the  notation  A{u)  for  the  design  equation  residual,  i.e., 

.4(u)  =  - Ll{<i>{u),u)X{u)  -  Fu(<?^(u),u)  (2.6) 

where  <l){u)  and  X{u)  in  (2.6)  are  solutions  of  the  state  and  costate  equations. 

2.3  The  Sensitivity  Gradients 

If  the  state  and  costate  equations  are  satisfied,  then  the  variation  of  the  cost  function  is 
given  by  (see  Eqn.(2.4)): 

SF=^-{u,A{u))^.  (2.7) 

This  equation  implies  that  the  gradient  of  the  functional  with  respect  to  the  control  vari¬ 
ables  is  given  by  —A{u): 

V,F{u)=--A{u).  (2.8) 

Therefore,  a  perturbation  of  the  control  variables  with  the  control  residuals  multiplied  by 

a  small  parameter,  namely  u  =  £>l(u),  will  result  in  a  reduction  of  the  cost  function  by 

SF=-e\\A\\l+0{e^).  (2.9) 


3 


2.4  Discretization 


When  discretizing  the  problem  it  is  possible  either  to  derive  the  necessary  conditions  for  a 
minimum  in  the  continuous  formulation  and  then  discretize  or  to  discretize  the  functional 
together  with  the  state  equation  and  then  derive  the  discrete  necessary  conditions.  In  the 
latter  ca.se  the  discrete  minimization  problem  is  given  by: 

min„h  on  r'‘ 

=  0  on 

As  the  grid  mesh  size,  /i,  goes  to  zero,  solutions  of  both  approaches  should  converge  to  the 
differential  solution.  However,  for  finite  mesh  size  discretization  and  necessary  conditions 
do  not  necessarily  commute.  The  solutions  of  both  should  be  within  the  discretization  error 
from  the  differential  solution.  For  simplicity  in  this  paper  we  used  the  first  possibility.  The 
discrete  state,  costate  and  design  equations  are: 

=  0  on 

+  =  0  on  O'*  (2.11) 

=  0  on  r\ 

We  define  A^{u'')  similarly  to  (2.6). 

2.5  A  Gradient  Descent  Algorithm 

The  following  is  a  gradient  descent  minimization  algorithm  which  follows  immediately  from 
the  above. 

1.  Start  with  an  initial  approximation  for  the  control,  Uq . 

2.  Solve  the  state  equation  for  <f)^ . 

3.  Solve  the  costate  equation  for  . 

4.  Compute  the  amplitude  of  the  perturbation,  ji ,  with  a  line  search, 

and  update  the  control  variables:  <—  u'^  +  0 

5.  If  the  residuals  of  the  state,  the  costate  and  the  control 
equations  are  greater  than  some  preassigned  value,  in  Lz  norm, 
then  goto  2;  else  stop. 

Note  that  steps  2,  3  and  5  consist  of  a  global  computation  over  the  whole  domain. 

The  complexity  of  this  algorithm  is  given  by  0{M^N‘),  where  M  is  the  number  of 
control  parameters,  N  is  the  number  of  grid  points,  and  p  and  /  are  integers  which  depend 
on  the  problem  and  the  PDE  solver  which  is  used  to  solve  the  state  and  restate  equations. 
For  example,  if  a  MG  solver  is  used  to  solve  the  PDFs  then  /  =  1. 
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3  A  Multigrid  One  Shot  Minimization  Method 

The  gradient  descent  algorithm  is  applied  on  a  sequence  of  nested  grids,  where  each  coarse 
grid  accelerates  the  convergence  rate  of  its  finer  grid.  On  each  grid  two  processes  are 
employed:  relaxation  and  coarse  grid  correction  of  all  the  variables,  including  the  control 
variables.  On  coarse  grids  the  state,  the  costate  and  the  design  equations  are  restricted 
from  the  finer  grid  with  the  full  approximation  scheme  [3]. 

3.1  Relaxation 

On  each  level  a  relaxation  is  performed  on  the  state,  costate  and  control  variables.  The 
state  and  costate  equations,  which  are  elliptic  PDEs,  are  relaxed  by  a  Gauss-Seidel  or 
damped  Jacobi  relaxations.  The  control  variables  are  relaxed  by 

^  u'*  +  (3.1) 

where  of  and  are  chosen  to  guarantee  good  smoothing  for  the  control  variables,  2is 
discussed  in  Sec. 4,  and  where  A'''{u^)  are  the  residuals  of  the  design  equation.  This  step 
should  be  followed  by  an  update  of  the  state  and  costate  solutions.  The  construction  of 
and  is  done  so  that  the  boundary  data  is  updated  with  a  high  frequency  dominated 
quantity. 

In  elliptic  systems  a  perturbation  of  the  boundary  condition  with  a  Fourier  mode 
has  an  exponential  decaying  effect  on  the  interior  solution  of  the  form  where  y  is  the 

distance  from  the  boundary  and  (t{u})  is  a  positive  monotonically  increasing  function  of  u> 
for  large  |aj|,  [1].  For  the  Laplace  equation  the  decaying  rate  is  given  by  Therefore,  in 

an  MG  scheme  it  is  preferable  to  perturb  the  boundary  condition  with  only  high  frequency 
modes  relative  to  the  given  level.  In  that  case  only  local  relaxations  will  be  needed  in 
order  to  update  the  solutions  after  each  optimization  step,  resulting  in  an  order  0{N~) 
operations  for  one  optimization  step.  N  is  the  number  of  interior  grid  points  on  a  given 
level,  and  d  is  the  space  dimension.  On  the  coarsest  grid  the  relaxation  of  the  control 
variables  is  given  in  Sec.  2.5  The  PDEs  are  solved  over  the  whole  domain  thus  taking  into 
account  the  lowest  frequencies.  In  that  way  the  set  of  grids  is  complete  in  the  sense  that 
all  Fourier  frequencies  are  treated  at  some  level. 

3.2  The  Coarse  Grid  Equations 

The  restriction  of  the  necessary  condition  for  a  minimum  to  the  coarse  grid  is  done  with 
the  full  approximation  scheme  (see  appendix). 

Coarse  Grid  State  Equation 

on  D" 


(3.2) 


where  and  are  restriction  operators  which  are  defined  on  the  interior  grid  points  of 
the  domain  and  which  are  not  necessary  identical  [.'Jl. 

Coarse  Grid  (Jostate  Equation 


LfX«  +  F«  =  /« 

A*)  =  if /«A*  +  /,f  A\  /«u*)  “  ^ 

-/,!'[i;'‘A'‘  +  Ff^'‘,A'‘,B‘y 

(Joarse  Crid  design  equation 


ifA"  +  F«  =  /« 

fH  =  I +  X\n'‘) 

A‘,u'‘)  =  if /«A‘  +  /"A*,  ;»!,'■) 

-/,f[i;‘A'‘  +  F„V.A'‘,«'')] 


on 


pH 


(3.4) 


where  and  are  restriction  operators  which  are  defined  on  the  boundary  r*‘,  and 
where  the  right  hand  sides  and  are  zero  on  the  finest  grid. 


3.3  The  Coarse  Grid  Cost  Function  and  Gradient 

It  can  be  shown  that  the  full  approximation  scheme  coarse  grid  equations,  (3. 2-3. 4),  are  the 
necessary  conditions  for  a  minimum  of  the  following  constrained  minimization  problem; 

min^H  F^{u^,<f>^)  -  {f^ ,<t>^)w  -  UH on  P" 

L»{<i>^,u»)  =  J^  onn», 

where  /^,  and  are  defined  in  Eqns.(3.2),  (3.3)  and  (3.4).  This  implies  that  the 
coarse  grid  gradient  is  given  by 

At" f  (if  A" +  F«). 

Thus  the  relaxation  defined  by  (3.1),  on  coarse  grids,  converges  to  the  solution  of  the  coarse 
grid  problem, 

X''(u'')  =  /.''-(i:''A«+F»).  (3.7) 

3.4  The  One  Shot  Minimization  Algorithm 

The  problem  is  solved  in  one  application  of  an  FMG  solver.  The  FMC  scheme  uses  a  Vcycle 
scheme  in  order  to  solve  the  problem  on  each  level.  The  Vcycle  is  composed  of  recursive 
applications  of  a  relaxation  and  coarse  grid  correction.  In  the  following  the  relaxation, 
Vcycle  and  FMG  schemes  are  presented. 
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Relaxation 

A  relaxation  sweep,  TV*,  is  defined  by  the  following: 

1.  Perform  one  relaxation  of  the  state  equation  for  (}>'\ 

2.  Perform  one  relaxation  of  the  costate  equation  for  A^‘ . 

3.  Update  the  control  variables  with  the  design  equation  residuals, 

uh  ^ 

4.  Perform  a  few  local  relaxations  of  the  state  and  costate 
equations  in  a  narrow  strip  near  the  boundary. 

Vcycle 

The  following  is  a  V^{ux^  V2)  cycle  where  V\  and  U2  are  integers  (in  the  numerical  demon¬ 
strations  we  used  i/j  =  2  and  1/3  =  O-  The  initial  grid  is  the  finest,  with  a  mesh  size  h. 

1.  Perform  i/j  relaxation  sweeps,  . 

2.  Restrict  the  state,  the  costate  and  the  design  equations 

to  the  coarse  grid  (Eqns . (3.2) , (3 .3)  and  (3.4),  with  H  =  2k) . 
Rescale  h  —*2h. 

3.  If  the  coarsest  level  is  not  reached  goto  1. 

4.  Solve  the  problem  with  the  standard  minimization  algorithm 
in  Sec.  2.5. 

5.  Interpolate  the  coarse  grid  correction  to  the  finer  grid: 


i — 

1 

+ 

on  ft*. 

A''  <- 

A'*  +  ^{A^'*  - 

on  ft*. 

on  r*. 

Rescale  | . 

Perform  1/3  rel2ucation  sweeps. 

7e*. 

7.  If  the  finest  grid  is  reached  then  stop,  else  goto  5. 

FMG  cycle 

The  following  is  a  n-FMG(t^i,  V2)  cycle  to  solve  the  problem  with  M  grids.  The  coarsest 
mesh  size  is  denoted  by  he- 

1.  Start  with  the  coarsest  grid,  (h  —  he) ,  and  solve  the  problem  with 
the  standard  minimization  algorithm  in  section  2.5. 

2.  Interpolate  the  solution  to  a  finer  grid,  rescale  h-~* 

3.  Perform  n  times  V^{v\,U2)  cycles. 

4.  if  the  finest  grid  is  reached  then  stop,  else  goto  2. 
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The  computational  cost  of  the  cycle  is  0{N)  operations,  and  it  reduces  the  error  of  the 
state,  costate  and  control  variables  by  an  order  of  magnitude  (in  L2  norm). 

4  Fourier  Analysis  of  the  Convergence  Rate 

Fourier  analysis  of  the  minimization  algorithm  is  described  next.  The  evolution  of  high 
frequency  errors,  in  the  control  variables,  is  considered  in  half  space.  Then  in  a  standard 
procedure  the  problem  in  half  space  is  reduced  to  the  boundary.  A  relation  between  errors 
and  residuals  of  the  design  equation,  on  the  boundary,  is  derived.  With  this  relation  the 
relaxation  and  coarse  grid  correction  of  the  control  variable  are  analyzed. 

4.1  Reduction  to  a  Boundary  Problem 

We  assume  that  the  state  and  costate  equations  are  satisfied  when  the  control  variables  are 
updated.  We  are  interested  in  the  amplification  factor  of  the  error  in  the  control  variables 
as  a  result  of  this  process.  In  the  vicinity  of  the  boundary,  the  non-smooth  errors  can  be 
analyzed  using  half  space  geometry.  This  approximation  is  valid  since  in  elliptic  problems 
non-smooth  Fourier  modes  decay  exponentially  into  the  interior  (see  Sec.  3.1).  Consider  a 
two-dimensional  geometry,  where  the  a:-axis  is  parallel  to  the  boundary  and  the  j/-axis  is  in 
the  normal  direction.  The  errors  of  the  state  and  costate  variables  satisfy  a  homogeneous 
equation  in  the  interior  at  every  optimization  step,  namely 

u^{x)  = 

where  (t{6)  and  <7{6)  are  determined  by  the  interior  state  equation: 

[^h^ixe/h^-a(9)y/h  _  Q 
j^h*^ix9/h^-^{e)y/h  _  Q 

By  substituting  these  expressions  into  the  boundary  conditions  of  the  state  and  costate 
error  equations,  we  obtain  relations  between  <^^(^,  j/  =  0),  X^{9,y  —  0)  and  u^(0),  which 
are  all  boundary  quantities.  Thus,  a  reduction  to  a  boundary  problem  has  been  obtained. 
From  the  boundary  problem  we  can  deduce  a  relation  between  the  residuals  of  the  design 
equation  and  the  errors  in  the  control  variables: 

A'^i9)  =  f'^{9)u'\9).  (4.3) 

T^{9)  is  the  symbol  of  the  Hessian  of  the  cost  function,  F,  subject  to  the  PDF  constraint. 
This  symbol  determines  the  smoothing  properties  of  the  control  variables  relaxation  as  well 
as  the  effectiveness  of  the  coarse  grid  correction,  as  is  discussed  next.  Note  that  the  explicit 
form  of  the  operator,  T^,  i^  not  known,  and  in  general  is  a  non-lccal  operator.  However, 
the  computation  of  its  symbol  is  straightforward. 


(4.1) 


(4.2) 
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4.2  The  Relaxation 


From  Eqns.(3.1)  and  (4.3)  it  follows  that  the  relation  between  the  errors  in  the  control 
variables  before  and  after  the  relaxation  is  given  by 

(4.4) 

where  the  relaxation  symbol  R^{6)  is  given  by 

=  (4  5) 

For  multigrid  proposes  it  is  desirable  for  R'*(d)  to  have  small  values  in  the  high-frequency 
range,  (f  <  |0|  <  tt).  In  that  case  the  relaxation  will  reduce  effectively  the  high-frequency 
errors  of  the  control  variables  prior  to  restricting  its  values  to  the  coarse  grid. 

Choice  of  High  Pass  Filter 

In  some  cases  the  relaxation  without  the  use  of  a  high  pass  filter  (HPF),  does 
not  smooth  the  errors  effectively  for  any  choice  of  0'^.  In  that  case  an  HPF  is  introduced 
as  a  preconditioner  of  the  control  residuals.  If  chosen  properly,  the  symbol  IF^{9)T'^{0)  is 
dominated  by  the  high  frequencies,  and  a  proper  choice  of  will  result  in  good  smoothing. 
The  HPF  is  particularly  effective  for  problems  in  which  the  transformation  T^{0)  is  a 
monotonically  decreasing  func  .ion  which  has  small  values  in  the  high  frequencies.  Without 
the  use  of  a  proper  HPF,  high-frequency  oscillatory  errors  might  enter  the  control  variables 
during  the  computation. 

Evaluation  of  the  optimization  step  size 

In  a  multigrid  cycle  the  relaxation  should  be  effective  mainly  in  the  high-frequency 
range.  The  relaxation  parameter  P'^  is  chosen  so  that  the  maximum  of  |.R^(0)|  in  the  high 
frequencies  will  be  minimal,  that  is, 

min  max  \\  +  p*^P^f’'{0)\.  (4.6) 


One  can  show  that  if  the  symbol  T'^(0}  does  not  change  sign,  then  p'‘  is  given  by 

2 


r  =  -73 


(^'‘T'*)m.n  +  (J^'*l'‘)r 


(4.7) 


where  and  {!F’^T^)max  are  the  minimal  and  maximal  values  of  T'^{0)T^{0)  range 

(f  <  1^!  <  ^)-  most  practical  problems  the  symbol  P^T^{0)  is  monotone.  Thus  p'^  is 
given  by 

o 

ah  _ _ 7 _  f4 


4.3  Two  Level  Analysis 

The  process  for  solving  the  optimal  control  problem  is  equivalent  to  a  process  of  solving  tlie 
equation  ,  =  r^,  where  eJ’’  and  are  the  errors  and  residuals  of  the  design  equation, 
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under  the  assumption  that  the  state  and  costate  equations  are  satisfied.  Using  standard 
multigrid  arguments  the  two  level  convergence  matrix  is  given  by  (See  Appendix) 

M'\e)  -  (4.9) 

where  R^{d)  is  the  relaxation  matrix 


R'^ie)  = 


1  -b  0 

0  1  +  +  tt) 


(4.10) 


and  C'‘(0)  is  the  coarse  grid  correction  matrix  given  by 


C\9)  = 


f''*(0  +  7r) 


The  asymptotic  convergence  rate  is  given  by  the  maximum  eigenvalue  of  the  matrix  M^{9), 
where  9  is  in  the  range  0  <  j0|  <  tt. 


5  Numerical  Examples 

In  this  section  we  demonstrate  the  performance  of  the  multigrid  one  shot  algorithm  and 
apply  the  analysis  developed  in  Sec. 4,  on  a  series  of  test  problems.  The  problems  are  solved 
in  a  two-dimensional  domain  which  is  defined  by 

n  =  {(x,y)|:0<a:<  1  ;  0  <  y  <  l}. 

The  constraint  is  the  Poisson  equation  and  the  boundary  conditions  are  periodic  in 
the  x-direction  and  Dirichlet  on  the  lower  boundary,  y  =  0.  The  minimization  problem  is 
defined  or  the  upper  boundary,  y  =  1. 

In  subsection  5.1  we  solve  an  optimal  control  problem  of  the  Dirichlet  boundary  con¬ 
dition  with  four  different  discretizations.  The  purpose  of  this  example  is  to  study  the 
dependence  of  the  convergence  behavior  on  the  choice  of  discretization.  Asymptotic  two 
level  convergence  rates  are  estimated  with  Fourier  analysis  for  the  different  discretization 
schemes.  The  predicted  and  the  actual  convergence  rates  are  compared. 

In  subsection  5.2  we  solve  an  optimal  control  of  the  Neumann  boundary  con  lition  for 
two  different  boundary  conditions.  The  purpose  of  this  example  is  to  show  the  use  of 
the  HPF  to  achieve  an  efficient  smoother  for  the  control  variables.  The  two  boundary 
conditions  correspond  to  qualitatively  distinct  transformations  between  error  and  residuals 
of  the  control  equation.  The  two  level  analysis  is  used  to  determine  a  proper  HPF. 

5.1  The  Dirichlet  Boundary  Control  Problem 

(Consider  the  minimization  problem  is  defined  by 

min  /  (^  ~ /*(^))  dx  +  T)  f  u^dx,  (5-1) 

«(x)  ^an  '  Jy=\ 
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where  rj  \s  a.  fixed  non  negative  parameter,  /*(x)  is  a  given  function  and  where  <p  satisfies 
the  state  equation 

A<p  —  f  on  n 

<  (f)  =  u{x)  on  y  =  1  ('">■2) 

(f>  =  (f>(i  on  y  =  0. 

It  is  easily  verified  that  the  the  costate  equation  is  given  by 

AA  =  0  on  il 

i  A  +  2(|J  - /‘(x))  =  0  on  y  =  l  (5,3) 

A  =  0  on  y  =  0 

and  the  design  equation  is  given  by 

.  dX 

•A  =  - - =  0  on  y  =  I.  (5.4) 

on 

5.1.1  Discretization 


We  have  used  four  different  discretizations  for  the  minimization  problem.  For  three  dis¬ 
cretizations  all  unknowns  were  defined  on  the  vertices  the  grid  lines  (referred  to  as  the 
“vertex  grid”).  The  control  variables  are  defined  on  the  intersections  of  the  grid  with  the 
boundary.  In  the  “cell  centered  grid”  the  variables  are  defined  on  the  centers  of  the  grid 
cells.  The  grid  is  extended  out  of  the  domain  and  virtual  cell-centered  points  are  defined 
on  the  neighboring  exterior  of  the  domain.  A  Dirichlet  boundary  condition  is  given  for  the 
average  of  the  variables  neighboring  the  boundary.  The  control  variables  are  defined  on  the 
centers  of  the  segments  connecting  the  intersection  of  the  grid  with  the  boundary.  Note 
that  in  the  multigrid  scheme,  the  vertices  of  the  grids  on  different  scales  are  nested  while 
in  the  cell-center  case  the  cells  faces  are  nested. 

In  the  vertex  grid  we  use  three  different  approximations  for  the  normal  derivative  on 
the  boundary: 

1)  A  first  order  approximation  for  the  normal  derivative 


d<j>  _  <f>i,2  —  <^1,1 
dm  h 


2)  A  second  order  approximation  for  the  normal  derivative 


(5.5) 


d4>  _  +  2^1,2  —  ^<^1,3 

dui  h 


(5.6) 


3)  A  use  of  a  virtual  point  out  of  the  domain,  were  its  Vd.lue  is  determined  with  the 
application  of  the  interior  operator  on  the  boundary 


1/  VI  “  *^’>1 

^  dn,~  2h 


A  cell  centered  discretization 

CC  : 


d4> 

dn. 


-  K-k 


—  *2 


(5.7) 


(5.8) 
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5.1.2  Reduction  to  the  Boundary 

In  the  following  we  analyze  the  design  equation  of  the  Dirichlet  boundary  control  problem 
in  the  discrete  space.  We  use  a  second  order  finite  difference  approximation  of  the  Lapiacian 
given  by 


(5.9) 


In  that  case  and,  therefore,  (t{9]  =■  o-(9).  The  term  in  Eqn.(4.1)  .satisfies  the 

following  second  order  equation  (see  Eqn.(4.2)) 

e''W  +  (_4-}-2cos0)4t"''‘'’'  =  U.  (.5.10) 

The  Fourier  Symbol  of  the  Normal  Derivatives 

The  normal  derivatives,  which  appear  in  the  design  equation,  have  the  following  Fourier 
symbols  for  the  different  discretizations: 


VXl  : 

II 

_  j 

h 

(■5 

11 

(  ^ 

V  y\.  Cl  * 

h 

VX  'i  : 

II 

g-o(5)  _  g<T(e) 

2h 

(5 

CC: 

II 

h 

(5 

The  Fourier  Symbol  of  the  Design  Equation 

In  terms  of  the  normal  derivatives  the  transformation  T  '{0)  (see  Eqn.{5.4))  is  given  by 


r‘(«)  =  -2[(£w)'  +  4  (S.I5) 

In  this  case  the  calculation  of  the  amplitude  of  the  minimization  step,  given  by  Eqn.(4.8) 
reduces  to 

/S'  =  - 5 - b - , - .  (.5,16) 

(£(fi) +(£w) +2'/ 

In  Fig.l  the  relaxation  symbol  k^{0)  =  1  -f  /T*T'^‘(0)  is  plotted  for  the  above  four 
discretizations.  For  all  four  discretizations  the  relaxation  reduces  the  high  frequency  errors 
by  a  factor,  smaller  than  0.5. 

Fig.2  depicts  the  maximal  eigenvalue,  |A|,„ox,  of  the  convergence  matrix  (4.9)  as  a 
function  of  the  number  of  minimization  steps,  i/,  on  a  given  level.  The  factor  by  which  the 
error  is  reduced  as  a  result  of  a  two  level  multigrid  cycle  is  bounded  by  |A(„,ax.  It  is  implied 
by  Fig.2  that  the  *  ell-centered  (CC)  and  second  order  vertex  (VX2)  schemes  are  expected 
to  have  a  better  performance  than  the  other  vertex  schemes. 
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Figure  1:  The  symbol  of  the  control  variable  relaxation  for  the  Dirichiet  boundary  control 
problem  with  r/  =  0. 


Figure  2;  Two  level  analysis  of  asymptotic  convergence  rates, |A|max,  as  a  function  of  the 
number  of  optimization  steps,  u,  for  t;  =  0. 


5.1.3  Convergence  Performance 

In  the  numerical  tests  the  problem  (5.1)-(5.2)  was  solved  for  the  four  discretizations  (5.5)- 
(5.8).  In  this  problem  there  was  no  need  to  use  a  high  pass  filter  since  tlje  transformation 
T^{9)  is  dominated  by  the  high  frequencies  in  all  four  discretizations.  The  minimization 
step  amplitude,  jS^.,  given  by  Eqn.{5.16)  was  used  in  the  computations.  The  rnuitigrid  one 
shot  algorithm  was  tested  using  between  two  and  seven  levels.  The  two  levels  convergence 
is  compared  with  the  convergence  predicted  by  the  analysis.  In  all  the  test.s  the  residuals 
of  the  state,  the  costaie  and  the  design  equations  were  computed  in  L-i  norm. 

In  the  two  levels  test  the  finest  grid  was  composed  of  2^  x  2^  grid  points  and  the  coarsest 
grid  was  composed  of  2®  x  2®  grid  points.  The  parameter  rj  was  set  to  zero.  In  Fig,3  the  two 
level  analysis  and  the  actual  convergence  rates  are  compared  and  the  similarity  between 
them  is  well  apparent. 

In  the  multilevel  test  the  fine  grid  was  composed  of  2"*  x  2”*  points,  with  m  =  5, 6,  7,  and 
the  coarsest  grid  was  composed  of  2  x  2  grid  points.  The  tests  with  different  choices  of  m 
were  done  in  order  to  check  if  the  algorithm  is  mesh  size  dependent.  All  the  results  in  Fig. 4 
were  done  with  a  cell-centered  discretization.  Since  the  case  =  0  in  (5.1)  corresponds 
to  a  trivial  problem,  the  case  rj  =  1  was  tested,  although  in  principle  the  results  should 
not  be  different.  Fig.4  A  shows  the  convergence  performance  of  the  analysis  problem  (5.3). 
Figs. 4  B  and  C  show  the  convergence  performance  of  the  optimization  problem  (5.1)  with 
rf  =  0  and  rj  —  I,  respectively.  The  depicted  residuals  in  4  B  and  C  are  the  average  of  the 
computed  state,  costate,  and  design  equations  residuals. 

In  all  problems  the  error  was  reduced  in  each  Vcycle  by  an  order  of  magnitude,  where 
each  Vcycle  has  a  computation  complexity  of  0{N). 
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VXITLA  O 
VXl  Numerical  + 


VX2TLA  O 
VX2  Numerical  + 


C.R.  0.25 


0  0.5  1  1.5  2  2.5  3  3.5  4  4.5 


0.5  1  1.5  2  2.5  3  3.5  4  4.3 


VX3TLA  O 
VX3  Numerical  + 


CCTLA  O 
CC  Numerical  -f 


C.R.  0.25 


0  0.5  1  1.5  2  2.5  3  3.5  4  4.5 

V 


0.5  1  1.5  2  2.5  3  3.5  4  4.5 


Figure  3:  Two  level  convergence  rates  (C.R.)  of  the  Dirichlet  boundary  control  problem 
as  a  function  of  minimization  steps  on  the  fine  level,  v.  “TLA”  stands  for  the  two  level 
analysis  prediction  and  “Numerical”  stands  for  actual  convergence  rate.  The  four  figures 
correspond  to  different  discretization  schemes  given  by  (5.1 1)-(5.14). 
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5.2  The  Neumann  Boundary  Control  Problem 

In  the  following  examples  we  compute  the  Neumann  boundary  condition  on  the  boundary 
y  =  1  such  that  the  values  of  the  solution  ^  on  that  boundary  will  match  some  given 
function. 

The  minimization  problem  is  defined  by: 

mm  f  U-f*{x)ydx,  (5.17) 

u(t)  Jy=l  ''  / 

where  (f>  is  satisfying  the  state  equation 

=  /  on  Q 

-  U  =  F(u)  on  y=l  (5.18) 

<f>  =  (po  on  3/  =  0. 


We  study  the  above  problem  for  two  different  right  hand  sides  of  the  boundary  condition 


Fi(«)  =  u 
F2{u)  =  Uj:. 


It  is  easily  verified  that  the  the  costate  equation  is 


(5.19) 


< 


AA  =  0 


A  =  0 


on  fl 
on  y  =  1 
on  y  =  0. 


(5.20) 


The  design  equation  is  given  on  the  boundary  y  =  I  for  the  two  right  hand  sides,  in  the 
corresponding  order,  by 


A  =  -A  =  0 
•^2  =  A^  =  0. 


(5.21) 


Discretization  The  state  and  costate  variables  were  discretized  on  a  cell  centered  grid. 
The  control  variables  are  defined  on  the  centers  of  the  segments  connecting  the  intersection 
of  the  grid  with  the  boundary  for  the  first  c<ise:  F^{u)  =  u.  In  the  second  case,  F-ziu)  =  u^, 
the  control  variables  are  defined  on  the  intersections  of  the  grid  lines  with  the  boundary. 


5.2.1  Analysis 

The  symbols  of  the  transformations  T^(^)  in  (4.5)  and  the  proper  amplitudes,  calculated 
with  Eqn.(4.8)  for  the  different  boundary  conditions  in  (5.19)  are  given  by 


'J'h(n\  _  _ 'n/^  2coa(S)  _  o  _  4 

'‘l  2  l-co«9  > 

f‘(«)  =  -V6-2cos«  ; 


(5.22) 
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One  can  immediately  observe  that  in  the  first  problem  the  transformation  symbol,  'fiid), 
is  not  dominated  by  the  high-frequencies.  This  means  that  no  choice  of  can  result  in 
a  relaxation  with  good  smoothing  properties,  since  high-frequency  errors  will  have  a  low 
weight  in  the  residuals  of  the  design  equation.  For  example,  ~  approaches  zero 

with  grid  refinement;  therefore  one  should  expect  very  slow  convergence  for  high-frequency 
errors.  For  this  case  the  two  level  analysis  predicts  a  non -con verging  scheme,  i.e.,  the 
maximal  eigenvalue,  |A|,„aa:,  of  the  convergence  matrix  (4.9)  is  greater  than  one. 

In  the  second  problem  there  is  a  high-frequency  dominance  in  the  symbol  T^iO).  There¬ 
fore  one  could  expect  that  a  few  local  relaxations  near  the  boundary  are  needed  after  each 
perturbation  of  the  control  variables.  In  this  case  the  two  level  analysis  predicts  the  same 
convergence  rate  as  in  the  analysis  problem,  see  Fig.6B. 

5.2.2  Use  of  a  High  Pass  Filter 

In  the  first  test  problem,  where  Fj{u)  =  u,  we  can  perturb  the  boundary  data  with, 
instead  of  A,  where  Dxx  is  a  second  order  tangential  derivative.  In  this  case  the  symbol  of 
the  design  equation  will  become 

DxMT’m  =  \/6-2cos{0),  (5.23) 

resulting  in  a  high  frequency  dominant  symbol.  The  relaxation  parameter  changes  and  is 
given  by 

A  use  of  a  higher  order  operator  such  as  a  fourth  order  tangential  derivative,  F  =  Dixxxi 
results  in 

Dxxxx{0)f!^{0)  =  -^v^6-2cos(0)(l  -  cos(^))  (5.25) 

with 

li>i{Dxx:xx)  =  2^.  (5.26) 

The  two  level  analysis  gives  a  much  better  convergence  performance  for  the  first  choice,  see 
Fig.6A. 

5.2.3  Convergence  Performance 

The  convergence  performance  of  the  Neumann  boundary  control  problem  is  tested  for  the 
two  cases  of  boundary  condition. 

In  Fig.7A,  the  convergence  of  the  optimal  control  problem  (5.17)-(5.18),  with  F{u)  =  u, 
is  depicted.  It  is  clear  from  Fig.7A  that  when  using  a  HPF  of  the  form  =  Du  the 
converegnce  rate  is  better  than  that  achieved  when  using  T  =  Dxxxx,  as  predicted  by 
the  analysis  (Fig.6A).  Without  using  a  HPF,  {F‘  =  /),  the  algorithm  didn’t  converge  and 
high-frequency  oscillatory  errors  were  observed  to  dominate  the  solution. 


A.  f'\{u)  =  u 


B.  =  u 


Figure  5;  The  symbol  of  the  control  variables  relaxation,  for  the  Neumann  boundary  control 
problem  (5,17)-(5.18).  In  Fig.A.  two  different  HPF  are  used,  ((5.23)  and  (5.25)),  with  their 
appropriate  relaxation  parameter  fS'‘,  (5.24)  and  (5.26). 

The  two  level  analysis  predicts  that  if  the  boundary  condition  is  changed  to  F{u)  —  Uj. 
in  Eqn.(5.18)  than  no  HPF  is  needed  for  the  problem  to  converge.  Fig.TB  shows  that  this 
is  indeed  the  case. 

6  Conclusion 

We  have  developed  a  multigrid  one  shot  algorithm  to  solve  the  infinite  dimensional  optimal 
control  problem.  An  analysis,  which  is  based  on  the  reduction  of  the  problem  to  Tie 
boundary,  was  performed  both  to  predict  the  convergence  rate  of  a  two  grid  algorithm  and  to 
determine  the  minimization  step.  Thus,  an  expensive  line  search  on  every  minimization  step 
was  not  required  on  fine  levels.  Numerical  demonstrations  on  a  series  of  two  dimensional 
test  problems  were  performed.  In  each  test  problem  the  amplitude  of  the  the  minimization 
step  on  fine  levels,  ,61^,  and  a  proper  high  pass  filter  (HPF),  was  determined  using 
the  analysis.  Comparison  of  the  two  level  convergence  rates  and  two  level  analysis  shows 
agreement  within  a  reasonable  error.  We  find  this  analysis  a  simple  and  powerful  tool.  In 
each  problem,  the  minimum  was  reached  at  a  cost  of  solving  the  analysis  problem  just  a 
few  times,  independent  of  grid  size. 
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Figure  6:  Two  level  analysis  convergence  rates  for  two  choices  of  boundary  conditions  in 
the  Neumann  boundary  control  problem  (see  Eqn.(5.19)).  T  stands  for  the  choice  of  high 
pass  filter. 


Figure  7:  A.  Convergence  rates  for  V(2, 1)  multigrid  cycle  of  the  Neumann  boundary  control 
problem,  with  two  different  HPFs:  Du  and  Dxxxx>  Fig.  B.  No  use  of  HPF  is  made.  Fine 
grid  contains  2^  x  2^  and  coarsest  4  grid  points. 
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A  A  Short  Review  of  Multigrid  Methods  for  Solving 
PDEs 

The  main  idea  of  multigrid  is  to  solve  the  problem  on  a  set  of  nested  grids.  On  each 
grid  there  are  two  processes:  relaxation  and  coarse  grid  correction  (CGC).  The  relaxation 
smooths  the  errors  while  the  CGC  eliminates  effectively  the  smooth  errors.  In  that  way 
a  very  efficient  algorithm  is  achieved  which  reaches  the  discretization  error  in  0{N)  op¬ 
erations,  where  N  is  the  number  of  interior  grid  points.  In  the  following  we  refer  to  the 
solution  of  a  discrete  elliptic  PDE  given  by 

iV  =  /\  (A.l) 

where  h  is  the  mesh  size. 

Coarse  Grid  Problem 

For  non-linear  problems  the  coarse  grid  equations  are  approximated  by  the  full  approx¬ 
imation  scheme  (FAS): 

=  /«/  +  t£(/),  (A.2) 

where  and  denote  the  restriction  and  interpolation  operators,  respectively,  and  where 
t[‘  is  defined  by 

(A.3) 

where  and  I{^  are  not  necessary  the  same  [3].  The  coarse  grid  correction  (CGC)  is  given 

by 

- /Mw)-  (A..) 

Smoothing  Properties 

In  the  infinite  space  mode  analysis  the  Fourier  symbol  of  the  differential  operator, 
L’\  is  denoted  by  For  standard  MG  to  work  properly  one  must  have  h-elliptic 

discretization,  i.e. 

|i‘(«)l  >  C|j|’"  /or  ||«|  <  ^  (A.5) 

where  L^{0)  is  the  Fourier  symbol  of  and  where  C  is  a  constant.  The  consequence  of 
the  condition  in  (A.5)  is  that  the  relaxation  will  be  effective  mainly  for  the  high  frequency 
errors.  For  the  Jacobi  relaxation,  the  relaxation  operator,  is  related  to  the  difference 
operator,  L^,  by 

R^  =  I-0L'\  (A.6) 

where  I  is  the  identity  operator  and  /?  is  a  parameter.  The  smoothing  rate  of  the  relaxation 
is  determined  by 

max  <  Co  <  1,  (A. 7) 

where  the  constant  Co  is  the  predicted  smoothing  factor  (typical  value  is  0.5). 

Two  Level  Analysis 
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For  the  linear  case  a  two  level  mode  analysis  can  he  done  to  predict  tin*  convergence 
rate  of  the  multigrid  cycle.  The  full  cycle  symbol  is  give  by 

M{9)  =  -  i\{e)L^ {29)-' 111 [9)V^[e)]k^(9Y\  (a.h) 

where  Rl’^  stands  for  the  -elaxation  operator,  I  is  the  unit  matrix,  stands  for  the  coarse 
grid  operator,  and  u\  and  Ui  are  integers. 

The  convergence  rate  of  the  cycle  is  give  by 

sup  \\M{9)1  (A.9) 

0<|(9|<f 
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